/****************************************************************************** 

  Copyright 2013 Scientific Computation Research Center, 
      Rensselaer Polytechnic Institute. All rights reserved.
  
  The LICENSE file included with this distribution describes the terms
  of the SCOREC Non-Commercial License this program is distributed under.
 
*******************************************************************************/
#include <PCU.h>
#include "crvShapeFixer.h"
#include "crvTables.h"
#include "maTables.h"
#include "maMesh.h" // for rotateTet
#include <pcu_util.h>

namespace crv {

/* projects vertex 3 onto the plane
   of the bottom triangle and returns
   the zone in which it lands as a bit code.
   Each bit indicates whether the area coordinate
   of that vertex is positive.
TODO: This needs to be updated to account for curved
elements. For now it is exactly like the one used in
linear mesh-adapt in ma/maShape.cc
*/

int getCrvSliverCode(
    Adapt* a,
    apf::MeshEntity* tet)
{
  ma::SizeField* sf = a->sizeField;
  apf::Mesh2* m = a->mesh;
  apf::Matrix3x3 J,Q;
  apf::MeshElement* me = apf::createMeshElement(m,tet);
  apf::Vector3 center(.25,.25,.25);
  apf::getJacobian(me,center,J);
  sf->getTransform(me,center,Q);
  J = J*Q; //Jacobian in metric space
  apf::destroyMeshElement(me);
  apf::Vector3 v03 = J[2];
  J[2] = apf::cross(J[0],J[1]); //face normal towards v[3]
  apf::Vector3	    projected = v03 - apf::project(v03,J[2]); //v[3] projected to face
  apf::Matrix3x3    inverseMap = apf::invert(apf::transpose(J));
  apf::Vector3	    basisPoint = inverseMap * projected;
  apf::Vector3	    areaPoint(1-basisPoint[0]-basisPoint[1],
			      basisPoint[0],
			      basisPoint[1]);
  int code = 0;
  for (int i=0; i < 3; ++i)
    if (areaPoint[i] > 0)
      code |= (1<<i);
  PCU_ALWAYS_ASSERT(code);
  return code;
}

ma::CodeMatch matchSliver(
    Adapt* a,
    apf::MeshEntity* tet)
{
  /* this table was auto-generated by the sliverCodeMatch program */
  ma::CodeMatch const table[8] =
  {{-1,-1}
  ,{4,1}
  ,{1,1}
  ,{2,0}
  ,{2,1}
  ,{0,0}
  ,{1,0}
  ,{0,1}
  };
  return table[getCrvSliverCode(a,tet)];
}

CrvFaceVertFixer::CrvFaceVertFixer(Adapt* a)
{
  mesh = a->mesh;
  edgeSwap = makeEdgeSwap(a);
  nes = nf = 0;
  edges[0] = 0;
  edges[1] = 0;
  edges[2] = 0;
}

CrvFaceVertFixer::~CrvFaceVertFixer()
{
  delete edgeSwap;
}

void CrvFaceVertFixer::setTet(apf::MeshEntity** v)
{
  /* in this template, the bottom face and v[3]
    are too close, the key edges are those that bound
    face v(0,1,2) */
  apf::findTriDown(mesh,v,edges);
}

bool CrvFaceVertFixer::requestLocality(apf::CavityOp* o)
{
  return o->requestLocality(edges, 3);
}

bool CrvFaceVertFixer::run()
{
  for (int i=0; i < 3; ++i)
    if (edgeSwap->run(edges[i]))
    {
      ++nes;
      return true;
    }
  ++nf;
  return false;
}

int CrvFaceVertFixer::getSuccessCount()
{
  return nes;
}

CrvEdgeEdgeFixer::CrvEdgeEdgeFixer(Adapt* a)
  : doubleSplitCollapse(a)
{
  mesh = a->mesh;
  adapter = a;
  edgeSwap = ma::makeEdgeSwap(a);
  nes = ndsc = nf = 0;
  sf = a->sizeField;
  edges[0] = 0;
  edges[1] = 0;
}

CrvEdgeEdgeFixer::~CrvEdgeEdgeFixer()
{
  delete edgeSwap;
}

void CrvEdgeEdgeFixer::setTet(apf::MeshEntity** v)
{
  /* in this template, the v[0]-v[2] amd v[1]-v[3]
    edges are too close. */
  apf::MeshEntity* ev[2];
  ev[0] = v[0]; ev[1] = v[2];
  edges[0] = findUpward(mesh, apf::Mesh::EDGE, ev);
  ev[0] = v[1]; ev[1] = v[3];
  edges[1] = findUpward(mesh, apf::Mesh::EDGE, ev);
}

bool CrvEdgeEdgeFixer::requestLocality(apf::CavityOp* o)
{
  return o->requestLocality(edges, 2);
}

bool CrvEdgeEdgeFixer::run()
{
  for (int i = 0; i < 2; i++) {
    if (edgeSwap->run(edges[i])) {
      ++nes;
      return true;
    }
  }
  if (doubleSplitCollapse.run(edges)) {
    ++ndsc;
    return true;
  }
  ++nf;
  return false;
}

int CrvEdgeEdgeFixer::getSuccessCount()
{
  return nes + ndsc;
}

CrvLargeAngleTetFixer::CrvLargeAngleTetFixer(Adapt* a):
  edgeEdgeFixer(a),
  faceVertFixer(a)

{
  adapter = a;
  mesh = a->mesh;
  tet = 0;
  fixer = 0;
}

CrvLargeAngleTetFixer::~CrvLargeAngleTetFixer()
{
}

int  CrvLargeAngleTetFixer::getTargetDimension()
{
  return 3;
}

bool CrvLargeAngleTetFixer::shouldApply(apf::MeshEntity* e)
{
  if (!ma::getFlag(adapter, e, ma::BAD_QUALITY))
    return false;
  tet = e;

  if ( ! getFlag(adapter,e,ma::BAD_QUALITY))
    return false;
  tet = e;
  ma::CodeMatch match = matchSliver(adapter,e);
  if (match.code_index==EDGE_EDGE)
  {
    fixer = &edgeEdgeFixer;
  }
  else
  {
    PCU_ALWAYS_ASSERT(match.code_index==FACE_VERT);
    fixer = &faceVertFixer;
  }
  apf::MeshEntity* v[4];
  mesh->getDownward(e,0,v);
  apf::MeshEntity* rv[4];
  ma::rotateTet(v,match.rotation,rv);
  fixer->setTet(rv);
  return true;
}

bool CrvLargeAngleTetFixer::requestLocality(apf::CavityOp* o)
{
  return fixer->requestLocality(o);
}

void CrvLargeAngleTetFixer::apply()
{
  if (! fixer->run())
    ma::clearFlag(adapter, tet, ma::BAD_QUALITY);
}

int CrvLargeAngleTetFixer::getSuccessCount()
{
  return fixer->getSuccessCount();
}

CrvLargeAngleTriFixer::CrvLargeAngleTriFixer(Adapt* a):
  remover(a)
{
  adapter = a;
  mesh = a->mesh;
  tri = 0;
}

CrvLargeAngleTriFixer::~CrvLargeAngleTriFixer()
{
}

int CrvLargeAngleTriFixer::getTargetDimension()
{
  return 2;
}

bool CrvLargeAngleTriFixer::shouldApply(apf::MeshEntity* e)
{
  tri = e;
  apf::Downward edges;
  mesh->getDownward(tri, 1, edges);
  remover.setEdge(edges[0]);
  return true;
}

bool CrvLargeAngleTriFixer::requestLocality(apf::CavityOp* o)
{
  return remover.requestLocality(o);
}

void CrvLargeAngleTriFixer::apply()
{
  ma::clearFlag(adapter, tri, ma::BAD_QUALITY);
  return;
}

int CrvLargeAngleTriFixer::getSuccessCount()
{
  return 0;
}

CrvShortEdgeFixer::CrvShortEdgeFixer(Adapt* a):
  remover(a)
{
  adapter = a;
  mesh = a->mesh;
  sizeField = a->sizeField;
  shortEdgeRatio = a->input->maximumEdgeRatio;
  nr = nf = 0;
  element = 0;
}

CrvShortEdgeFixer::~CrvShortEdgeFixer()
{
}

int CrvShortEdgeFixer::getTargetDimension()
{
  return mesh->getDimension();
}

bool CrvShortEdgeFixer::shouldApply(apf::MeshEntity* e)
{
  if ( ! getFlag(adapter,e,ma::BAD_QUALITY))
    return false;
  element = e;
  apf::Downward edges;
  int n = mesh->getDownward(element,1,edges);
  double l[6] = {};
  for (int i=0; i < n; ++i)
    l[i] = sizeField->measure(edges[i]);
  double maxLength;
  double minLength;
  apf::MeshEntity* shortEdge;
  maxLength = minLength = l[0];
  shortEdge = edges[0];
  for (int i=1; i < n; ++i)
  {
    if (l[i] > maxLength) maxLength = l[i];
    if (l[i] < minLength)
    {
      minLength = l[i];
      shortEdge = edges[i];
    }
  }
  if ((maxLength/minLength) < shortEdgeRatio)
  {
    ma::clearFlag(adapter,element,ma::BAD_QUALITY);
    return false;
  }
  remover.setEdge(shortEdge);
  return true;
}

bool CrvShortEdgeFixer::requestLocality(apf::CavityOp* o)
{
  return remover.requestLocality(o);
}

void CrvShortEdgeFixer::apply()
{
  if (remover.run())
    ++nr;
  else
  {
    ++nf;
    ma::clearFlag(adapter,element,ma::BAD_QUALITY);
  }
}

}
